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Abstract. - The effect of polydispersity on the phase behaviour of hard spheres is examined 
using a moment projection method. It is found that the Boublik-Mansoori-Carnahan-Starling- 
Leland equation of state shows a spinodal instability for a bimodal distribution if the large 
spheres are sufficiently polydisperse, and if there is sufficient disparity in mean size between the 
small and large spheres. The spinodal instability direction points to the appearance of a very 
dense phase of large spheres. 



The phase behaviour of hard sphere mixtures has received much attention recently following 
the observation by Biben and Hansen p| that, with an improved closure approximation, inte- 
gral equation theory indicates the presence of a fluid-fluid spinodal instability in a bidisperse 
mixture if the ratio of the large to small sphere diameters is > 5. Several experimental groups 
report confirmatory evidence H |l 0]. The dense phase of the large spheres is sometimes 
observed to be an ordered phase and at other times to be an amorphous colloidal glass, 
although some discrepancies in detail remain |3[ |J. 

Numerical solutions to integral equations frequently only suggest the existence of a spinodal 
instability, from the apparent divergence of the structure factor for instance. Thus for binary 
hard spheres, Caccamo and Pellicane J^] can argue that the integral equation evidence supports 
fluid-solid segregation rather than fluid-fluid separation. Recent computer simulations support 
this scenario, by establishing that a region of fluid-fluid coexistence exists, but is metastable 
with respect to fluid-solid coexistence [pi. 

Analytic solutions to integral equations, although rare and generally perceived to be less 
accurate than numerical solutions, do not suffer from the same problems. Most notable of these 
for hard spheres is the Percus-Yevick closure 0. As reported by Lebowitz and Rowlinson [|| 
though, this closure predicts that a binary hard sphere fluid is stable for all size ratios and 
physically accessible compositions. 

The appeal of the Percus-Yevick approximation is tempered by an internal inconsistency 
between the compressibility equation of state (EOS) and the virial EOS. It is precisely this 
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inconsistency that the new closure approximations mentioned above are designed to cure. A 
phcnomcnological way around this problem though is to interpolate between the two EOS's. 
For monodispcrse hard spheres, a suitable interpolation leads to the highly succesful Carnahan 
and Starling EOS ||. For an arbitrary mixture, the same interpolation leads to an EOS first 
considered by Boublik and Mansoori et al (BMCSL EOS) (uj. This improved EOS also 
predicts phase stability for a binary hard sphere fluid for any size ratio and composition. 

It is well known though that polydispersity enhances phase instability. An interesting and 
experimentally relevant question therefore concerns the effects of size polydispersity on the 
phase behaviour of hard sphere mixtures. Fortunately the above EOS's have been extended 
to polydisperse systems, for instance the Helmholtz free energy and related thermodynamic 
quantities corresponding to the BMCSL EOS are given by Salacuse and Stell jnj in terms of 
the number density and the first three moments of the size distribution. 

For the Percus-Yevick compressibility EOS, Vrij |1| established that an arbitrarily polydis- 
perse hard sphere fluid is stable, thus confirming a conjecture of Lebowitz and Rowlinson 
However the question of whether, under the BMCSL EOS, a binary mixture becomes unstable 
if one or both of the species is allowed to become polydisperse has not to my knowledge been 
previously investigated. Perhaps surprisingly in the light of Vrij's result the answer to this 
question is "yes". Moreover the instability reflects the experimental observations, albeit at 
a much larger degree of polydispersity and size ratio than are seen in reality. Very recently 
Cuesta has found a similar spinodal instability in hard spheres with a unimodal, log-normal 
size distribution Jl3[ . 

I address the spinodal stability problem using a novel moment projection method developed 
recently by Sollich and Cates Q and myself Q . The approach reduces the free energy to one 
which depends only on moment densities (defined below). It gives exact results for cloud and 
shadow curves (the polydisperse analogues of the binodal) and for spinodal curves, provided 
the excess free energy is a function only of the corresponding moments of the size distribution. 
Since the BMCSL excess free energy is precisely of this form, the new tool may be aplied 
to this situation. By way of contrast, Cuesta approaches the spinodal problem by casting it 
as an integral equation |l3). He also finds that if the excess free energy depends only on a 
few moments, the integral equation is reducible to a matrix problem, which can in fact be 
proved to be identical to the results obtained below. Such a reduction in dimensionality for 
the spinodal problem has been noted several times in the past [fl6f . 

I outline the approach in general first to indicate its application to spinodal curves, which 
was not previously covered in detail. Suppose that the excess free energy depends only on a 
few moment densities, defined to be <f> n = fni&^/V, where V is the system volume and 

the functions / n (c) are arbitrary. The sum is over all particles in a homogeneous phase, and 
the &i are individual particle properties, such as size. If we set /o(c) = 1 then we include the 
number density p = 4>q amongst the moment densities (in doing this though the n = term 
should be omitted from certain sums below). As previously reported, moment densities can 
be treated as independent thermodynamic density variables provided the ideal part of the free 
energy density is suitably generalised: 

/< id ) = pk B T (In p-l)-pTs(m n ) (I) 

where s(m n ) is a generalised entropy of mixing per particle, and m n — <fi n /p — X)j=i fn{&i)/N 
(n > 0) are generalised moments. An expression for s(m n ) can be given only if the parent or 
feedstock distribution, p(<r), is specified. This is because one must know the parent distribution 
from which particles are drawn, to conclude anything about the distributions in daughter 
phases. The generalised entropy of mixing is then given by a Legendre transform (I will set 
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dh 

s(m n ) = h(9 n ) + 9 n m n , m„ = - — , (2) 

n>0 



fca = T = 1 in the following for simplicity): 

$ T"n Tr) — — 

where the function h is a generalised cumulant generating function for p(cr), 

h = ln / dap(a) exp|^- ^ 9 n f n (a) . (3) 

« ~, ^ n 



n>0 



The excess free energy density, f^\4> n ), is added to /( ld ) to obtain the overall free energy 
density /. The spinodal stability condition is then the standard one of positive definiteness of 
the matrix of second partial derivatives of / with respect to the </>„ (including (j)Q — p amongst 
these). For the ideal part these second derivatives are readily shown to be 



1 1 d 2 s d 2 f ld) lv- 9 2 s 

V "i„m r - , — = - > m r - , (n>0), 

n p — ' om n om r ooaocpn p ^-^ om n om r 

n,r>0 r>0 

9 2 /( id) 1 Q2 S 



a 2 / (id) i i Q2s d 2 f^_ 

>n 

d(j> n d(f> r pdm n dm r ' ' ^ ^ ^ 

The key is clearly the matrix d 2 s/dm n dm r It can be easily shown from eq. (^|) that if 
(A)„ r = d 2 h/d9 n d9 r (n, r > 0), then d 2 s/dm n dm r = — (A _1 )„ r . 

Since the parent distribution enters into the solution in an essential way, it constrains the 
mean values of the moments to those of the parent distribution. This constraint is essential to 
obtain the exact spinodal and related stability conditions, and cloud and shadow curves. For 
a homogeneous system, it means that after all calculations have been completed the results 
should be constrained to lie on the physical dilution line, 4> n /p = m n = (/n(c)) p where 
(. . .) p — J da p(a)(. . .). Fortunately this constraint has a simple implementation. By explicit 
differentiation of eq. (|3j) the moments are given as functions of 9 n by 



dap(a)f n (a)exp\-^29 r f r (a) . (5) 

"* r>0 



Since the map 9 n <-> m n should be 1 : 1, otherwise the Legendre transform is ill defined, 
the dilution line constraint is seen to correspond to the point 9 n — 0. This observation now 
makes it obvious why stability conditions such as the spinodal condition, the critical point, 
and so on depend only on a finite set of moments or cumulants of the parent distribution. 
These conditions involve second or higher order derivatives of h evaluated on the dilution 
line, which corresponds to the point 9 n = 0. Since h is a cumulant generating function, 
these derivatives are precisely the cumulants of p{cr). This observation generalises a number of 



moment truncation theorems reported previously by other workers 16 . Applying this dilution 
line constraint to the matrix A obtains for instance 

(A) nr = (fn(<y)fr(<y))p - (fn(cr)) p (f r (a)) p . (6) 

This matrix is readily inverted for use in eqs. (Q). 

The spinodal stability limit corresponds to a vanishing eigenvalue of the matrix of second 
partial derivatives of the free energy. The associated eigenvector contains valuable but oft- 
neglected information on the spinodal instability direction. This is the direction in composition 
space in which the homogeneous system starts to become unstable towards small perturbations, 
as the spinodal line is crossed. It can be used to determine (mean-field) critical points, as 
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points where the instability direction is tangent to the spinodal line. Away from such points, 
the instability direction tends to reflect the slope of ticlincs connecting coexisting phases. 

In the Sollich-Cates picture [0, the spinodal instability is contained in the family of 
size distributions p(a)exp[^2 \ n f n (a)]. Expanding this about a point on the dilution line 
indicates that, in terms of a size distribution, the spinodal instability direction is contained 
in Ap(cr) = pp(a)^2 AX n f n (a), where the prefactor p comes from the n = term. The 
AA„ can be found from the constraints J da f n (a)Ap(a) = A(j> n where A(f> n is the prescribed 
eigenvector. From this, if (B)„ r = {f n (<y)f r (a)} p (n,r > 0), the spinodal instability direction 
is 



Actually, in the Sollich-Cates picture one can also show that the matrix <9 2 /( ld ) j 'dcf) n d(f) r — 
((pB) _1 ) nr [0. The equivalence to eqs. (^) can be established with some algebra. Note that 
if the /n(c) (n > 0) form an orthonormal set with p(cr) as a weight function, the matrices A 
and B become unit matrices, and many of these results simplify markedly. 

I now apply this machinery to examine the spinodal stability of hard spheres. It is convenient 
to introduce a fiducial diameter ctq and scale all densities with 7t<7q/6. The excess free energy 
density corresponding to the BMCSL EOS is given by Salacuse and Stell Jll| 



The moment densities are 4> n — {^crQ/6)pm n and the moment functions themselves are 
/n(c) = (cr/co)™- Thus (f>o = it<Jqp/6 is now a dimensionless number density, and (f>3 = 4> 
is the volume fraction. 

For the parent distribution I take either a Schulz distribution with an extra monodisperse 
component or a mixture of two Schulz distributions. I will characterise these by the mean 
diameter and the degree of polydispersity, where the latter is the ratio of the standard deviation 
to the mean diameter (usually expressed as a percentage). The moments ((a/ao) n ) p which 
appear in the matrices A and B are given by simple algebraic relations once the distributions 
have been specified. 

The 4x4 spinodal determinant is constructed as described above, and examined as a function 
of the ratio of mean diameters, cr^/as, the degree of polydispersity of the two components, 5l 
and 5s, and the volume fractions of the two components, (f>L and <j>s- Generally the determinant 
is always positive except when a large value of the diameter ratio, ul/^s si 60, is combined 
with a large degree of polydispersity of the large spheres, Sl ^ 50%. 

Fig. 1(a) shows the spinodal instability appearing as the diameter ratio is increased, for 
a bimodal distribution with a monodisperse small component. Fig. 1(b) shows the same 
for fixed mean diameter ratio and increasing degree of polydispersity of both components. 
Interestingly, at fixed mean diameter ratio and degree of polydispersity of the large component, 
the instability region diminishes as the degree of polydispersity of the small species increases. 
This is illustrated in fig. 2. In both figures, results are only shown for the physically reasonable 
region (f> s + 4>l < 0.5. 

The tick marks (short lines) on the spinodal lines indicate the spinodal instability direction 
projected into the (</>s,</>l) plane. The insets to these figures show unnormalised parent 
distributions at one selected point in each plot together with the spinodal instability direction 
from eq. (^). Thus in fig. 1(b) for example, at the marked point, the system is unstable 
towards a composition fluctuation in which the density of large spheres increases and the 
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density of small spheres decreases: the precise functional form of the composition fluctuation 
is the dashed line in the inset. 

The tick marks and insets show that the instability is one of demixing with some additional 
size partitioning. The instability direction, as indicated by the tick marks for instance, is never 
found to be tangential to the spinodal line in the physical region although it does approach this 
condition for high values of the large sphere volume fraction. This clearly points to a critical 
point lying at large 4>l and suggests that the equilibrium state is one of a highly dense fluid 
of large spheres coexisting with much less dense fluid of mostly small spheres. Note that the 
degree of polydispersity for the large spheres is easily large enough to suppress the formation 
of an ordered phase jl8) . 

Why should polydispersity have this effect? If one introduces small spheres into a fluid of 
large spheres, depletion forces in the latter are apparently insufficient to bring about phase 
separation, at least according to the BMCSL EOS. It is known though that polydispersity 
decreases the pressure of a hard sphere fluid |To| , due to relaxed packing constraints. It appears 
from the present calculation that polydispersity can make the large sphere fluid sufficiently 
compressible for the small sphere depletion forces to induce fluid-fluid phase separation. What 
is remarkable is that the BMCSL EOS appears to capture this effect, without being specifically 
designed to do so. 

These results show intriguing parallels with the experimental observations, with regard to 
the nature of phase separation in hard sphere mixtures. In the experiments ^ || 0], phase 
separation is observed at a much smaller diameter ratio, typically ul/<Js ~ 10, and degree 
of polydispersity, typically < 10%. The dense phase of large spheres is also observed to 
be crystalline or amorphous. These differences are presumably due to the sensitivity of the 
phase boundary to the accuracy of the EOS, as reported by Biben and Hansen Jj| (given this 
sensitivity, one can speculate whether the discrepancies between experimental results [|[ [| 
for size ratio dependence may be due to polydispersity effects). Also, if the polydispersity in 
the dense phase is sufficiently small, ordering into a crystal phase is likely to take place. The 
novelty of the present analysis is that the broad picture is reproduced by the BMCSL EOS, if 
polydispersity is taken into account. 

Some intriguing possibilities still remain though for the true equilibrium phase diagram of 
an arbitrarily polydisperse mixture. One such is fractionation into multiple crystalline phases 
[ fl9| . Another is to relax the packing constraints in a dense crystal by degrading the crystal 
symmetry from FCC, or the formation of superlattice structures which could be regarded as 
interpenetrating multiple crystals. At high volume fractions though, the glass transition and 
other kinetic considerations undoubtably play a major role in determining what is observed 
experimentally. 

I acknowledge useful correspondence with Jose Cuesta, Richard Sear and Peter 
Sollich. 
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Fig. 1. - Spinodal instability lines for polydisperse hard sphere mixtures (solid lines with tick marks 
in main plots): (a) monodisperse small spheres mixed with polydisperse large spheres, Sl = 50%, for 
fi/cs = 60, 80 and 100, and (b) mixture of polydisperse small and large spheres at ctl/o\s = 100, 
for Sl = Ss =50%, 65% and 90%. The tick marks indicate the spinodal instability direction. The 
insets show the (unnormalised) size distribution (solid lines) and the direction in which it is spinodally 
unstable (dotted lines), at the marked points. 
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Fig. 2. - Mixture of polydisperse small and large spheres at ul/us = 100, Sl = 50%, for Ss = 
(monodisperse), 30% and 60%. 



